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Abstract 

This  research  was  directed  towards  developing  new  knowledge  of  combustion  processes  in  liquid- 
propellant  rocket  engines.  Attention  was  given  to  the  combustion  of  different  fuels,  beginning  with  hydrogen- 
oxygen  but  ultimately  including  multicomponent  hydrocarbon  mixtures  representative  of  RP  fuels,  and  to 
nonlinear  processes  of  importance  in  combustion  instability.  Effects  of  detailed  chemical  kinetics  on 
performance  and  stability  were  addressed  by  introducing  systematically  reduced  chemistry  that  enables 
theoretical  analyses  to  be  completed  thoroughly  and  accurately.  The  interactions  between  chemistry  and 
turbulence  were  addressed,  including  the  effects  of  the  heat  release  on  the  turbulence.  The  orientation  of  the 
work  was  fundamental  and  aimed  at  predicting  not  only  performance  and  acoustic  response  but  also  intrinsic 
instability.  The  results  were  intended  to  help  to  improve  understanding  of  combustion  mechanisms  and 
combustion  instabilities  in  liquid-propellant  rocket  motors.  The  most  recent  specific  results  concern  induction 
times  in  hydrogen-oxygen  mixtures  and  heal-release  effects  on  turbulent  mixing. 

Research  Highlights 

Above  the  second  explosion  limit,  the  ignition  of  hydrogen  mixtures  is  known  to  exhibit  a  relatively 
long  period  of  chain-carrier  growth  prior  to  significant  heat  release,  which  occurs  through  chain  recombination 
once  a  sufficiently  large  radical  pool  has  formed.  Since  this  induction  period  typically  represents  a  considerable 
fraction  of  the  overall  combustion  time,  designs  of  hydrogen  combustion  devices  must  account  for  this  delay. 
Because  of  its  practical  importance,  the  induction  time  in  hydrogen-air  ignition  has  remained  an  active  research 
topic  for  over  40  years.  Numerous  shock-tube  experimental  measurements  of  hydrogen  induction  times  are 
available  [l]-[9].  The  criteria  used  in  experiments  for  defining  the  induction  time  vary.  Thus,  for  some  authors 
the  induction  period  ends  at  the  onset  of  luminosity  or  as  the  pressure  rises  above  a  given  threshold  value  as  a 
result  of  heat  release,  while  others  prefer  to  relate  the  induction  time,  t„  to  the  radical  concentration  through  OH 
emission  or  absorption.  Under  most  conditions,  in  contrast  to  ignition  of  many  other  fuels,  the  resulting  value  of 
t,  for  hydrogen  depends  only  weakly  on  the  definition  used. 

In  addition  to  the  experimental  work,  there  have  been  a  number  of  analytical  studies  devoted  to  the 
calculation  of  the  hydrogen-oxygen  induction  time  at  high  temperature  [1,  5,  6,  7,  10,  11].  Although  the  success 
of  many  of  these  previous  analytical  and  computational  efforts  has  been  limited  by  uncertainties  in  reaction-rate 
kinetics,  these  works  provided  much  understanding  on  the  induction  process.  For  instance,  Schott  and  Kinsey 
[1]  already  anticipated  that  the  radicals  O  and  OH  follow  a  steady-state  approximation  if  the  mixture  is 
sufficiently  rich,  a  result  that  Trevino  [11]  obtained  later  in  his  systematic  reduction  of  the  ignition  chemistry. 
The  induction  time  was  determined  in  [1]  from  the  evolution  of  the  H-atom  concentration,  the  only  species  out 
of  steady  state.  Both  reactant  consumption  and  heat  release  were  seen  to  be  negligible  during  this  initial  stage  of 
radical  growth,  giving  an  H-atom  concentration  that  increases  exponentially  in  time.  Since  H  +  O2  -4  OH  +  O  is 
the  rate-limiting  reaction,  the  resulting  induction  time  was  seen  to  be  proportional  to  the  reciprocal  of  the  initial 
oxygen  concentration  [1].  The  calculation  procedure  can  be  extended  by  accounting  for  a  larger  number  of 
radicals,  as  done  in  [5,  6,  7,  10],  to  give  a  system  of  linear  differential  equations  for  the  radical-pool  evolution. 
The  resulting  induction  time  becomes  inversely  proportional  to  the  largest  eigenvalue  of  the  associated 
eigenvalue  problem,  showing  a  dependence  on  the  mixture  composition  more  complex  than  that  of  rich 
mixtures. 

The  present  work  combines  the  knowledge  currently  available  on  the  H2-O2  kinetics  [12]-[15]  with  the 
procedure  outlined  above  for  the  calculation  of  t,  [5,  6,  7,  10],  thereby  providing  an  improved  theoretical 
prediction  for  the  induction  time,  valid  under  variable  conditions  of  pressure,  temperature  and  composition.  The 
kinetics  of  hydrogen  combustion  is  now  well  established,  and  only  small  differences  (mainly  in  reaction-rate 
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constants)  exist  between  the  different  mechanisms  currently  available  (e.g.,  [12]-[15]).  Use  of  these  chemical- 
kinetic  mechanisms  enable  the  calculation  of  ignition  histories  by  numerical  integration  of  the  conservation 
equations  for  species  and  energy  in  a  homogeneous  mixture  at  constant  pressure.  Results  of  such  detailed- 
chemistry  computations  are  compared  in  Fig.  1  with  the  experimental  measurements  of  Just  and  Schmalz  [7] 
and  Bhaskaran  et  al.  [8],  the  induction  time  being  identified  in  the  numerical  computations  with  the  peak  H- 
atom  concentration.  The  results  obtained  with  the  different  mechanisms  are  very  similar;  significant 
discrepancies  are  only  found  as  the  crossover  temperature  of  the  second  explosion  limit  [16]  is  approached. 
From  the  comparisons  with  the  experimental  results,  it  can  be  inferred  that  the  mechanism  of  Miller  and 
Bowman  [13]  tends  to  underpredict  the  crossover  temperature,  while  the  GRI  mechanism  [14]  and  that  compiled 
by  Maas  and  Wamatz  [12]  tend  to  overpredict  it.  For  the  four  conditions  of  pressure  and  composition  used  in 
the  figure,  the  San  Diego  mechanism  [15],  developed  in  the  present  project,  is  the  one  that  shows  the  best 
agreement  with  the  experiments  over  the  whole  range  of  temperatures  considered.  In  view  of  these  results,  the 
analytical  development  uses  the  reaction  rates  given  in  internet-based  updated  version  of  [15]. 

The  reaction  mechanisms  [12]-[15]  comprise  about  20  reversible  reactions  among  8  reactive  species, 
including  as  intermediates  H,  O,  OH,  HOj  and  H2O2.  The  numerical  computations  indicate  that  the  six  steps 
suffice  to  describe  accurately  the  induction  period.  The  rate  constants  vary  with  the  temperature  T  according  to 
ki  =AiT"‘  exp(-7;^/r)  [15],  where  the  units  are  mol/cm^  s'*  and  K.  Nonunity  chaperon  efficiencies  of  0.3  for 

O2  and  7.0  for  H2O  must  be  considered  when  calculating  the  third-body  concentration  Cm  appearing  in  the  rate  of 
reaction  5  [15].  Note  that  the  present  short  mechanism  does  not  consider  HO2  consumption,  a  simplification  that 


mav  introduce  small  inaccuracies  in  induction  times  very  close  to  crossover  [4]. 

No. 

Reactions 

A 

n 

T 

0 

H2  +  O2  OH  +  OH 

1.70  X  10*^ 

0.0 

24044 

1 

H2  +  02^  HO2  +  H 

2.91  X  10'^ 

0.47 

26800 

2 

H  +  02->  OH  +  0 

3.52  X  10*® 

-0.7 

8535 

3 

0  +  H2  ^  OH  +  H 

5.06  X  lO" 

2.7 

3145 

4 

OH  +  H2  ^  H2O  +  H 

1.17  X  10’ 

1.3 

1813 

5 

H-i-02  +  M^H02  +  M 

2.6  X  10*’ 

-1.2 

0 

As  can  be  seen,  two  different  initiation  reactions  are  considered  in  the  mechanism,  while  others  that  are 
much  slower,  such  as  H2  —>  2H,  are  neglected.  Except  for  an  early  investigation  [10],  all  previous  efforts  have 
used  only  one  of  the  two  reactions,  i.e.,  reaction  0  is  selected  in  [1,  5,  6]  while  reaction  1  is  preferred  in  [7,  11]. 
The  rate  constant  for  H2  +  O2  OH  -l-  OH  was  first  estimated  by  Semenov  [17]  and  later  improved  by  Ripley 
and  Gardiner  [18].  The  reaction-rate  parameters  given  above  for  0  are  those  proposed  by  Jachimowski  and 
Houghton  [19],  which  are  adopted  by  the  four  currently  used  mechanisms  [12]-[15].  The  rate  of  reaction  1  is 
determined  from  measurements  of  the  inverse  reaction  HO2  -H  H  — >  I^  -t  Q.  Experimental  data  collected  in  the 
eighties  [20,  21,  22]  have  been  used  to  determine  the  reaction-rate  parameters  used  for  this  reaction  in  [12]-[15]. 
The  different  constants  differ  only  moderately  over  the  range  of  temperatures  of  interest  (1000  K  <  T  <  4000 
K),  e.g.,  by  less  than  ten  percent  for  the  two  internet-based  mechanisms  [14,  15].  With  the  rate  constants 
proposed  in  [15]  reaction  1  is  the  dominant  chain-initiation  reaction  at  temperatures  below  T  “  2315  K,  and  the 
opposite  behavior  is  observed  above  this  temperamre.  Therefore,  both  reactions  must  in  principle  be  retained  in 
the  mechanism  for  increased  accuracy. 

The  analytical  development  addresses  three  differential  equations  for  the  time  development,  which  can 
be  written  nondimensionally  as 

V^)yH  3^0  )^OH  ^0’ 

d>'o/dT  =  (^/'CoXj'h  /2  -  ^o)  (2) 

and  ~  (^/^oh)(3’h/2  +  iVo  ""  J’oH 

The  asymptotic  analysis  results  in  explicit  simple  formulas  for  the  induction  times  for  large  and  small  values  of 
the  equivalence  ratio  <t>  and  a  more  complex  exact  solution  that  applies  for  all  0.  Figure  2  compares  these 
results.  It  is  seen  that  the  asymptotic  result  for  lean  mixtures,  which  has  H  in  steady  state,  is  good  for  ^  <  0.07 , 
the  asymptotic  result  for  stoichiometric  and  rich  mixtures,  which  has  O  and  OH  in  steady  state,  is  good  for 
0  >  1,  and  in  the  intermediate  range  0.07  <  0  <  1  all  three  radicals  are  out  of  steady  state,  and  the  more  complex 
exact  solution  is  needed.  Nevertheless,  the  sum  of  the  two  asymptotic  results  provides  induction  times  that  are 
correct  within  and  order  of  magnitude  over  the  entire  range  of  equivalence  ratio. 
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In  other  recent  work,  evaluations  were  made  of  DNS  (direct  numerical  simulation)  results  for  the 
mixture-fraction  scalar  Z  and  the  scalar  dissipation  %  in  temporally  evolving  shear  layers.  Effects  of  high  heat 
release  have  not  been  available  previously.  For  high  heat  release,  typical  of  hydrocarbon  combustion,  the 
mixing  was  found  to  be  substantially  different  than  without  heat  release.  The  probability-density  function  of  the 
scalar  and  the  conditional  rate  of  scalar  dissipation  were  found  to  be  affected  by  the  heat  release  in  such  a  way 
that  the  heat  release  substantially  decreases  the  overall  reaction  rate.  To  help  clarify  implications  of  the 
assumptions  underlying  popular  models  for  interaction  between  turbulence  and  chemistry,  the  local  structure  of 
the  scalar  dissipation  rate  at  the  reaction  sheet  was  extracted  from  the  DNS  database.  The  applicability  of 
flamelet  models  for  the  rate  of  scalar  dissipation  was  examined.  To  assist  in  modeling,  a  characteristic  length 
scale  was  defined,  representing  the  distance  around  the  reaction  sheet  over  which  the  scalar  field  is  locally 
linear,  and  statistical  properties  of  this  length  scale  were  investigated.  This  length  scale  can  be  used  to 
investigate  values  of  the  rate  of  scalar  dissipation  that  mark  the  boundary  between  flames  that  feel  a  constant 
scalar  dissipation  field  and  those  that  do  not. 

The  topology  of  isoscalar  surfaces  is  central  to  validity  of  assumptions  underlying  modeling  of 
nonpremixed  turbulent  combustion.  This  is  especially  true  for  flamelet  approaches  that  rely  on  expansions 
about  the  stoichiometric  value  Z  =  but  it  also  applies  to  other  approaches  that  deal  with  averages  of  % 
conditioned  on  Z.  The  simplest  conceivable  configuration  is  that  of  a  planar  stagnant  molecular  mixing  layer,  in 
which  isoscalar  surfaces  are  parallel  planes,  and  is  constant  everywhere,  where  D  is  the  scalar  diffusion 

coefficient.  Since  this  uniformity  of  |VZ|  from  Z  =  0  to  Z  =  1  is  unrealistic  in  turbulence,  steady  counterflow 
laminar  flamelet  models  have  been  considered,  which  maintain  parallel  planar  isoscalar  surfaces  but  have  |VZ| 
varying  smoothly  from  zero  at  Z  =  0  or  at  Z  =  1  to  a  maximum  near  Z  =  0.5  in  each  instantaneous  realization. 
Turbulence,  however,  generates  much  more  complex  topology  than  this  for  isoscalar  surfaces,  causing  them  to 
be  nonplanar  and  time-dependent,  a  fact  that  is  responsible  for  the  knovra  limitations  of  steady  laminar  flamelet 
models,  including  those  that  take  into  account  the  linear  variation  of  %  with  Z  near  Z  =  Z^.  Turbulence  can 
generate  more  than  mild  wrinkles  in  the  curves  on  space-time-resolved  bases.  Ridges  can  develop,  with  |VZ| 
achieving  minima  along  certain  surfaces,  as  well  as  zero-gradient  points,  where  The  wrinkles  readily 

produce  maxima  and  minima  in  instantaneous  plots  of  %  as  a  function  of  Z  along  paths  in  the  direction  of  VZ, 
and  the  minima  descend  to  =  0  at  zero-gradient  points.  Such  plots  must  in  fact  end  when  X-^’  turning  back 
on  themselves  in  at  least  one  other  maximum-gradient  direction,  at  maximum  or  minimum  values  of  Z  different 
from  1  or  0,  unless  the  zero-gradient  point  has  an  inflection  type  of  behavior  instead  of  being  a  local  absolute 
extremum.  These  complexities  deserve  consideration  in  modeling. 


A  first  step  in  this  direction  is  to  study  in  greater  detail  the  local  behaviors  of  the  Z  and  x  fields  in  the 
vicinity  of  isoscalar  surfaces.  Let  y  be  the  distance  in  the  direction  normal  to  the  surface  Z  =  Z^,  with  y  >  0  for 
Z  >  Z^  and  y  <  0  for  Z  <  Z^.  A  Taylor  expansion  through  second  order  in  y  then  is 


Z  =  Z^  +  a,y-l-fly, 

1  .  .A 

(4) 

where 

II 

(5) 

and 

1  (9'Z 
a,  =  —n,n.  ^  ^  , 

2  dx^dx^ 

(6) 

in  which 

1  dZ 

a,  aXi 

(7) 

is  the  i'th  component  of  a  unit  vector  normal  to  the  surface.  All  partial  derivatives  here  are  evaluated  the  point 
where  Z  =  Z^.  The  second  coefficient,  differs  from  the  curvature  of  the  isoscalar  surface  Z  =  Z^  itself, 
which  is 


dn,^\(  d^Z  d^Z  ' 
(Jc,  ^  dxf  '  ‘  dx^dxj  ^ 


(8) 
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Both  and  Ui  are  random  variables.  Since  ;^;  =  2Daf ,  there  have,  in  effect,  been  numerous  statistical 
studies  of  a^,  but  has  been  relatively  neglected.  The  expansion  of  x  about  Z  =  Z^,  however,  depends  on  02- 


Through  first  order  in  (Z  -  Z^) , 


dZ 


(z-z.). 


where 


dZ 


z=z, 


dy^  dZ 


dy. 


dD 


=  SDfl,  "I"  2 — a, , 


(9) 


(10) 


the  last  term  of  which  may  be  neglected  because  D  typically  reaches  a  maximum  very  near  Z  =  Z,.  The  good 
approximation 

1  dx\ 


(11) 

8D,  tiZ|z=z, 

thus  relates  to  the  slope  of  the  curve  of  x  as  a  function  of  Z  at  Z  =  Z^.  Since,  according  to  this  equation, 
stagnant  mixing-layer  flamelet  models  have  a^=0,  while  steady  counterflow  laminar  flamelet  models  with 
Z^  <  0.5  have  >  0,  the  statistics  of  pertains  to  the  potential  applicability  of  these  models. 

Another  reason  for  interest  in  is  its  relevance  to  limitations  on  the  use  of  the  characteristic  length 

f  =  Z,V2D7i;  =  Z>,  (12) 

often  considered  in  turbulence  modeling.  This  length  defines  the  characteristic  distance  over  which  Z  varies  by 
a  fractional  amount  of  order  unity  only  when  the  quadratic  term  is  small  compared  with  the  linear  term.  This 

condition  breaks  down  when  |y|  becomes  of  order 

=  (13) 


The  departure  of  Z  from  Z^  at  this  limit  is 

AZ  =  Z^d/i  =  al/\a,\,  (14) 

which  may  be  called  the  limiting  scalar  deviation  in  that  Z  ceases  to  vary  linearly  with  distance  y  beyond  this 
deviation.  It  thus  is  a  measure  of  the  amplitude  of  wrinkles  of  curves  of  Z  as  a  function  of  distance  normal  to 
isoscalar  surfaces.  This  limiting  deviation  and  the  associated  limiting  length,  d,  which  is  the  characteristic  linear 
dimension  or  length  of  such  wrinkles,  are  random  variables  having  statistics  determined  by  those  of  Oi  and  Ui- 
The  statistics  of  Ci  and  02,  conditional  on  Z  =  Z^,  thus  have  bearings  on  a  number  of  questions. 

Graphs  of  the  conditional  averages  of  %  as  functions  of  Z  were  investigated  in  this  work.  These 
averages  are  basically  equivalent  to  conditional  averages  of  Oi-  Corresponding  graphs  of  conditional  averages  of 
Ui  also  were  investigated.  For  heat  release  Q  =  0  they  are  antisymmetric  about  Z  =  0.5,  as  they  must  be.  They 
are  also  very  nearly  antisymmetric  for  Q  ^  0,  suggesting  that,  although  the  dynamics  correspond  to  Z^  =  0.2  in 
these  cases,  contrary  to  ai  the  asymmetry  of  the  density  field  has  little  effect  on  the  conditioned  mean  02,  thereby 
allowing  any  value  of  Z  to  be  interpreted  as  Z^.  The  density  change  has  little  effect  on  the  conditional  mean 
(fljlz)  for  0.3 <Z <0.7,  but  there  are  noticeable  influences  in  the  extremes  representative  of  hydrocarbon 
combustion.  The  conditioned  means  all  peak  near  Z  =  0.1,  but  the  peak  for  Q  =  0  is  much  higher  than  those  for 
Q  0.  The  heat  release  of  combustion  therefore  reduces  Cj ,  which  would  tend  to  increase  the  limiting  scalar 
deviation  and  the  wrinkle  length. 

It  is  found  from  the  probability-density  function  for  that  the  mode,  while  positive,  is  nearly  at  02  =  0, 
and  there  is  an  appreciable  probability  that  is  negative  in  all  cases.  This  means  that  (dx/dZ)^^^^  is  negative 


nearly  as  often  as  it  is  positive.  Instantaneous  laminar-flamelet  profile  shapes  of  the  scalar  dissipation  thus  are 
the  exception  rather  than  the  rule.  Instantaneous  profiles  of  x{^)  Reynolds  numbers  therefore  are 

highly  irregular,  with  (dx/dZ)^^^  =  0  likely  being  a  better  approximation  than  any  counterflow  flamelet  model. 
More  research  therefore  is  needed  in  modeling  x{^)  in  turbulent  combustion. 
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Figure  1:  The  lines  represent  the  variation  with 
temperature  of  the  induction  time  as  obtained  from 
detailed-chemistry  numerical  computations  (solid  hnes: 
SD  mechanism  [15],  dotted  lines:  Maas  and  Wamatz 
[12],  dashed  lines:  GRI  [14],  dot-dashed  lines:  Miller 
and  Bowman  [13]),  while  the  symbols  denote  the 
shock-tube  experimental  results  of  Just  and  Schmalz  [7] 
(p  =  0.41,  0.43,  1.41  atm)  and  of  Bhaskaran  et  al.  [8]  (p 
=  2.5  atm). 


Figure  2:  The  induction  time  at  atmospheric 
pressure  as  obtained  from  integrations  of  the 
full  conservation  equations  (solid  lines),  from 
evaluations  of  the  exact  solution  to  (l)-(3) 
(dashed  lines),  and  from  evaluations  of  the 
asymptotic  predictions  (dotted  lines). 
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Relevance  /  Transitions 

•  New  high-performance  rockets  need  liquid  propellants  for  specific  impulse,  control  and  throttling. 

•  Codes  predicting  performance  and  instability  can  give  wrong  results  because  of  wrong  chemistry. 

•  Correct  chemistry  is  too  complicated  for  codes. 

•  Reduced  chemistry  retains  essentials  of  correct  chemistry  but  is  simple  enough  to  use  in  codes. 

•  Result  is  improved  prediction  of  motor  performance  and  instability. 

•  Recent  applications  have  been  to  predicting  emissions  of  oxides  of  nitrogen  in  turbulent  combustors  using 

KTVA  code  (Solar  Turbines)  and  predicting  autoignition  of  natural  gas  in  dual-fuel  diesel  engines  (BKM 

Co.). 

•  Rocket  applications  pertain  to  space  launch,  long-range  missiles  and  high-altitude  maneuverability. 

•  Propellant  combinations  include  LOX-hydrogen  and  LOX-RP. 

•  Other  major  contributions  have  been: 

1 .  Discovery  of  bimodal  probability-density  function  for  combustion  instability  as  a  result  of  nonlinear 
subcritical  bifurcation,  explaining  why  liquid-propellant  rockets  can  sometimes  jump  from  stable 
burning  to  noisy  high-amplitude  acoustic  instability. 

2.  Explanation  of  Hewitt  empirical  correlation  for  stable  and  unstable  LOX/RP-1  rockets  with  like-on-like 
injector,  in  a  diagram  of  frequency  and  ratio  of  injector  diameter  to  injection  velocity,  on  the  basis  of 
acoustic  amplification  by  mixing-layer  combustion. 

3.  Identification  of  cause  of  droplet  burning  time  being  minimum  at  critical  pressure,  on  the  basis  of  three- 
component  phase  diagrams,  explaining  new  mechanism  of  combustion  instability  at  near-critical 
conditions,  with  more  stable  operation  at  subcritical  or  supercritical  conditions. 

4.  Development  of  two-step  to  four-step  reduced  chemistry  for  hydrogen-oxygen  and  hydrocarbon-oxygen 
rocket  combustion,  providing  accurate,  fundamentally  based  and  practically  useable  description  for 
calculating  liquid-propellant  rocket  performance  and  instability. 
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